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Abstract 

Using lattice simulations, we study the extent of the conformal window for an SU(3) gauge theory 
with Nf Dirac fermions in the fundamental representation. We extend our recently reported work, 
describing the general framework and the lattice simulations in more detail. We find that the 
theory is conformal in the infrared for Nf = 12, governed by an infrared fixed point, whereas the 
Nf = 8 theory exhibits confinement and chiral symmetry breaking. We therefore conclude that 
the low end of the conformal window Nj lies in the range 8 < Nf < 12. We discuss open questions 
and the potential relevance of the present work to physics beyond the standard model. 

PACS numbers: ll.10.Hi, 11.15.Ha, 11.25.Hf, 12.60.Nz 



I. INTRODUCTION 



The conformal window in a gauge field theory with Nf light fermions is the range of Nf 
values such that the theory is asymptotically free and the infrared coupling is governed by an 
infrared fixed point. In an SU(iV) gauge theory with Nf Dirac fermions in the fundamental 
representation, the conformal window extends from UN/ 2 down to some critical value N'j at 
which a transition is expected to a phase in which chiral symmetry is broken spontaneously, 
and confinement sets in. In a recent paper [I], we provided the first nonperturbative evidence, 
using lattice simulations, that the lower end of the conformal window for the SU(3) gauge 
theory lies in the range 8 < Nj < 12. 

Gauge theories in or near the conformal window could play a key role in describing physics 
beyond the standard model. If the fermions carry electroweak quantum numbers, and if 
Nf lies outside but near the conformal window, then the theory could drive electroweak 
breaking, forming the basis of walking technicolor theories. If the fermions do not carry 
electroweak quantum numbers, then Nf could lie within the conformal window, and the 
theory could describe some new, conformal sector, possibly coupled to the standard model 
through SU(iV) invariant operators. It is important to learn as much as possible about the 
extent of the conformal window in these theories, as well as the order of the transition at 
Nf and the properties of the theory within the window and near it. 

To obtain the result 8 < Nj < 12 for Dirac fermions in the fundamental representation of 
an SU(3) gauge group, we employed pQ a gauge invariant, nonperturbative running coupling 
derived from the Schrodinger functional of the gauge theory j2j [31 Hj. Defined within a 
Euclidean box of volume 0(L 4 ), it avoids typical finite volume effects by using L itself as 
the sliding scale. For the asymptotically free theories being considered, it agrees with the 
perturbative running coupling coupling at small enough L, and can be used to probe for 
conformal behavior in the large L limit. We made use of staggered fermions as in Ref. [3], 
and therefore restricted attention to values of Nf that are multiples of 4. The value Nf = 16 
leads to an infrared fixed point that is so weak that it is best studied in perturbation 
theory. The value Nf = 4 is expected to be well outside the conformal window, leading to 
confinement and chiral symmetry breaking [6J as with Nf = 2. We thus focused on the two 
values Nf = 8 and Nf = 12. We argued |T] that for Nf = 12, the theory is indeed conformal 
in the infrared. For Nf = 8, we showed, in disagreement with an earlier lattice study [7], 
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that the theory breaks chiral symmetry and confines. There is no evidence for an infrared 
fixed point. 

In this paper, we provide a more detailed description of the results of Ref . [T] , and extend 
the analysis in several ways. Continuing to focus on the values Nf = 8 and 12, we describe 
more extensive numerical simulations of the running coupling and discuss in more detail the 
treatment of lattice artifacts and the extrapolation to the continuum. We again conclude 
that 8 < < 12, with a more precise determination of the large L behavior of the running 
coupling. 

This work paves the way for future SU(3) simulations at other values of Nf, in particular, 
in the range between 8 and 12, for fermions in other representations of the SU(3) gauge group, 
and for other gauge groups. Simulations using other definitions of the running coupling, for 
example, derived from the Wilson loop, should also be carried out [8]. Conclusions about 
the conformal window based on the study of running couplings should be confirmed by zero- 
temperature lattice simulations of the chiral condensate and other quantities. These include 
the particle masses and Goldstone boson decay constants for Nf just below NS, and various 
correlation functions within the conformal window [9]. 

In Sec. [TTJ we describe what is known from perturbative and other studies about the con- 
formal window in SU(iV) gauge theories. For comparison, we describe briefly the conformal 



window in supersymmetric SU(iV) gauge theories. In Sec. Ill, we review the Schrodinger 



functional framework [21 El S] for our numerical simulations. Our lattice simulations are de- 



scribed in Sec. [TVJ We discuss the algorithms, the use of staggered fermions, the continuum 
extrapolation, and analysis methods. In Sec. [V] we present our results for both Nf = 8 and 
12, and compare to other studies. We summarize our work and discuss open questions in 



Sec. VI The details of our error analysis are given in Appendix [Aj and tables of simulation 



data appear in Appendix [Bj 



II. THE CONFORMAL WINDOW 



We first review what is known from the perturbative expansion of the beta function, and 
then discuss briefly a nonperturbative approach based on the counting of degrees of freedom 
and some other, quasiperturbative studies. 
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A. Perturbation theory 



The existence of a conformal window in SU(iV) gauge theories has been known since the 
computation of the two-loop beta function by Caswell in 1974 [TU]. If the number of massless 
fermions Nf is near but just below the number Nf f at which asymptotic freedom sets in, 
then the two-loop term is opposite in sign to the one-loop term and the resultant infrared 
fixed point is weak, accessible in perturbation theory. There is no confinement and chiral 
symmetry is unbroken. The properties of this phase were studied by expanding in Nf { — Nf 
in Ref. [II]. 

As Nf is reduced, the strength of the infrared fixed point grows, with Nf ultimately 
reaching the value Nj at which the transition to the chirally broken and confining phase is 
thought to set in. There is no a priori reason to expect the infrared fixed point to remain 
perturbative through this window, although arguments to this effect have been advanced 

If the theory is formulated in the continuum and a running coupling g 2 {L) is defined at 
some length scale L, we have L(d/dL)g 2 (L) = (3 (g 2 (L)), where 

P (g 2 (L)) = ht{L) + b 2 f{L) + b 3 f(L) + b,g w (L) + • • • . (1) 

For the case of SU(3), the first two (scheme-independent) coefficients are 

bi = j£f2 [11 - §#/] M = j^i [102 - f N f ] • (2) 

The next two coefficients depend on the renormalization scheme. In the MS scheme, they 
are given by [13J 
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bf s = jj— (29243.0 - 6946.30A/ + 405.089A) + 1.49931AJ) . (4) 
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and 

2 

For Nf close to 33/2, the two-loop infrared fixed point value g 2 is very small, and therefore 
corrected very little by the higher order terms. 

If the loop expansion is reliable to estimate gl, other quantities can also be estimated. An 
example is the (scheme-dependent) parameter 7 governing the approach to the fixed point. 
If the beta function is linearized in the vicinity of the fixed point, 

-p(f{L))^ 1 [gl-f{L)] , (5) 
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then as L — > oo, the approach to the fixed point from either side is given by 

f(L)^t- C -^. (6) 

For Nf = 12, there is a two-loop infrared fixed point at g 2 ~ 9.48, corrected to ~ 5.47 
at three loops in the MS scheme, and to ~ 5.91 at four loops. The critical exponent is 
then 7 = 0.36 at two loops, and in the MS scheme is given by 7 = 0.296 at three loops 
and 7 = 0.282 at four loops. The convergence of the loop expansion is not guaranteed, but 
the fact that the expansion parameter at the fixed point gl/iir is relatively small suggests 
that it could be reliable, and therefore that Nf = 12 lies inside the conformal window. For 
Nf = 8, there is no two-loop infrared fixed point. A fixed point can appear at three loops 
and beyond in some schemes, but its scheme dependence and typically large value means 
that there is no reliable evidence for an infrared fixed point accessible in perturbation theory. 
A nonperturbative study is essential. 



B. An upper bound on Nf 

We next review a conjectured inequality which leads to an upper bound on Nf [14]. For 
any asymptotically free theory, the thermodynamic free energy may be computed perturba- 
tively as T — > 00, approaching the (free) Stefan-Boltzman expression — (7r 2 T 4 /90)/uv, with 
Ajv = Wb + (7/8)4 A F ], where N B and N F are the numbers of (underlying) bosonic fields 
and four- component Dirac fields. Similarly, as T — > 0, if the effective low energy theory is 
infrared free, the free energy approaches the expression — (7t 2 T 4 /90)/ir, where /ir counts 
the (massless) infrared degrees of freedom in the same way. The conjectured inequality is 
simply / IR < /uv- 

For a nonsupersymmetric SU(A) theory with Nf massless Dirac fermions in the fun- 
damental representation, /uv — 2(iV 2 — 1) + (7/8)4NNf. If this theory is in the chirally 
broken phase at zero temperature, then /ir simply counts the number of Goldstone bosons: 
/ IR = Nf — 1. The inequality then gives N c f < A^/N 2 - 16/81. For SU(3), this is consistent 
with Nf = 12 being within the conformal window. 

It is interesting to note that for a supersymmetric SU(A r ) gauge theory with Nf massless 
Dirac fermions in the fundamental representation, where Nf denotes the transition point 
between the phase with infrared conformal symmetry and the free-magnetic phase, the 
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same inequality gives Nj < (3/2)iV, a limit precisely saturated by the result from duality 
arguments [IS]- It is natural to ask to what extent the inequality is also saturated in the 
nonsupersymmetric case. 

C. Other studies 

Finally, we note that several groups over the years have attempted to determine the value 
of Np as well as the nature of the transition as Nf Nj from below and features of the 
bound-state spectrum in this limit, by studying continuum Schwinger-Dyson (SD) equations 
[T6| ITT1 fl8] . Some truncation of the SD equations must be adopted. It is then assumed that 
the infrared behavior is governed by an infrared fixed point appearing in the perturbative 
beta function, and solutions corresponding to broken chiral symmetry are sought. This 
leads to a value for Nj slightly below 4N, approaching it in the large-iV limit, as well as 
information about the theory in the broken phase near the transition. The reliability of 
these results is not clear, however, since higher order effects are not obviously small. 

III. THE SCHRODINGER FUNCTIONAL FORMALISM 
A. Introduction 

The Schrodinger functional is the partition function describing the quantum mechanical 
evolution of a system from a prescribed state at time t = to another state at time t = T 
in a spatial box of size L with periodic boundary conditions [21 El E] • Dirichlet boundary 
conditions are imposed at t — and t = T where T is O(L). They are chosen such that 
the minimum-action configuration is a constant chromo-electric background field of strength 
0(1/ L). This can be implemented in the continuum |2J or with lattice regularization [19]. 
In either case, by considering the response of the system to small changes in the background 
field, a gauge invariant running coupling can be defined, valid for any coupling strength. 

The Schrodinger functional can be represented as the path integral 



(7) 
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where A is the gauge field and ip, ip are the fermion fields. W and W are the (fixed) 
boundary values of the gauge fields, and (, (, ( are the boundary values of the fermion 
fields at t = and t = T, respectively. 

Although the Schrodinger functional can be formulated completely in the continuum, 
from here on we will introduce a Euclidean spacetime lattice. The quantity Sg is chosen to 
be the standard Wilson gauge action [20] with a boundary improvement counterterm: 



p 



^ £ TtF^ - A(i _ Ct y Yfi{t - 0) + 8(t - T)]Tr F 0u F^ + 0{a" 



where Tr represents a color trace, a is the lattice spacing, and (3 = 2N/g$ with g the 
lattice coupling constant. The improvement coefficient w(P) = c t when multiplying timelike 
plaquettes which intersect the Dirichlet boundaries, and is equal to 1 elsewhere. For this 
computation we set c t equal to its value as determined in one-loop lattice perturbation theory 
0, 

ct = 1 + ^[-0.08900(5) + 0.00474(l)iV7]. (9) 
The operator F M „ is defined similarly to the continuum field strength tensor, 

F^p = A^A,, - AlA^ + [An, A v ) (10) 

with A* g{x) = (g(x + afi) — g(x))/a the discrete forward derivative operator. If we take the 
continuum limit of the action ([8]), we recover the standard Yang-Mills action. The sum over 
plaquettes P may be expanded out in terms of individual gauge links: 

J^Re Tr U P = ^ J^Re Tr [U u (x)U^x + u)Ul(x + ft + i>)[/J(a: + p)] . (11) 

For the fermionic action, we use the staggered approach as in Ref. [5], which reduces 
the 16 doubler species of a naively discretized fermion field to 4 degrees of freedom. In the 
continuum limit, a single staggered fermion field can be interpreted as four degenerate Dirac 
fermion fields. For Nf divisible by four, the total fermionic action Sf is then given by 

s F = Yl 5 /> ( 12 ) 

N f /i 
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where Sf is the fermion action for a single staggered field as in [5], 

Sf = \ Yl ^( x )x( x ) [ u n( x )x(x + A) - Ul(x - ji)x(x - //)] , (13) 

with 77^ the usual staggered phase factor rj^ = (— iy-' v <i* Xv . 

Without affecting the action, the spatial-periodicity condition can be generalized to in- 
clude a phase rotation on the fermion fields at the spatial boundaries, 

X (x + Lk) = e ie * X (x), x( x + Lk) = x(x)e~ ie \ (14) 

where k runs over all of the spatial directions. Imposing a nonzero value on the Q k has been 
shown in QCD to improve the ratio of the largest and smallest eigenvalues of the Dirac 
matrix [21], which in turn improves computational speed. However, this result is based on 
nonperturbative tuning, and is not guaranteed to extend to the theories we are considering. 
We therefore set 6^ = for simplicity. 



B. Temporal boundary values and definition of the Schrodinger functional cou- 
pling 

The gauge boundary values W, W are chosen such that the minimum-action configuration 
is a constant chromoelectric field whose magnitude is of 0(1/ L) and is controlled by a 
dimensionless parameter rj. The Schrodinger functional (SF) running coupling is then defined 
in terms of the response of the action to variations in rj. The setup is as follows: we take for 
the particular boundary values of the gauge fields 



U k (x, t = 0) = exp(aCfc) 
U k (x, t — T) = exp(aC' k ) 

where the Ck,C' k are spatially constant and abelian, 



(15) 
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(16) 



Classically, boundary conditions of this form lead to a constant chromoelectric back- 
ground field, with field strength proportional to 1/L. We adopt the particular set of bound- 
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1, 



— 7] — 71 



(17) 



ary values 

^3 = — 2^ + 3 " '", - 2' I A 

which are chosen to ensure that the background field is a stable solution to the classical field 
equations for small rj [22] . 

With the boundary conditions fixed in this way, the SF running coupling g 2 (L,T) is 
defined by taking 

(18) 



9 (L, T) 
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where 



k = 12 
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The factor k is chosen so that g 2 (L,T) equals the bare coupling at tree level. In general, 
g 2 (L, T) can be thought of as the response of the system to small changes in the background 
chromo-electric field. 

The fermionic Dirichlet boundary values (, £, (' are subject only to multiplicative renor- 
malization for staggered fermions [23J. As we are free to choose these values, we take them 
equal to zero, simplifying the calculation. 

The staggered approach to discretization of fermions can be thought of as splitting the 
16 degrees of freedom of a single spinor over a 2 4 hypercube of lattice sites. This framework 
makes it evident that staggered simulations require an even number of lattice sites in each 
direction. Thus with periodic boundary conditions in space, the spatial extent L/a of 
the lattice must be even. However, in the Schrodinger functional formalism, the Dirichlet 
boundaries in the time direction require an odd temporal extent T/a in order for the number 
of lattice sites to be even, since the sites located at t = and t = T are distinct. 

As a result, one cannot simulate with T = L, only with T = L±a. In the continuum limit 
T = L is recovered, but at a finite lattice spacing this results in the introduction of 0(a) 
lattice artifacts into observables. This is particularly undesirable, since staggered fermion 
simulations contain bulk artifacts at 0(a 2 ) and higher. Fortunately, simulating at T = L±a 
and averaging over the observed values eliminates this effect 0. We adopt this technique 
here, defining the central observable 



g\L) 
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y (A L — a) g 2 (L,L + a 



(20) 



which depends on only one large distance scale L. To be more explicit, this running coupling 
can be written as g 2 (/3,L/a) where f3 = 2N/g 2 . From this point on we will fix iV = 3, and 
so (3 = Q/g 2 . 



C. Schrodinger functional perturbation theory 

The SF coupling g 2 (L) has been normalized to give the bare lattice coupling g\ at tree 
level. With the lattice as a regulator, it can be expanded as a power series in g 2 with 
coefficients depending on a/L. By rearranging this series in terms of a coupling defined at 
an arbitrary scale and setting to zero terms that vanish as a — > 0, a continuum beta function 
can be defined. Its perturbation expansion leads to the universal coefficients b\ and b 2 of 
Eq. g at the one- and two-loop levels. 

The three-loop, scheme-dependent coefficient has been computed in the Schrodinger func- 
tional scheme by combining the two-loop perturbative computation of the SF running cou- 
pling in lattice perturbation theory with a similar lattice computation of the MS coupling 
constant [19]. The result is 



lSF iMS , ^2 C 2 fol( c 3 - cj) ( -,\ 

^ ~ ^ + 47T 16^ ' W 



where bf s is given by Eq. g with c 2 = 1.256+0.040^/ and c 3 = c^+1.197(10)+0.140(6)A^ / - 
0.0330(2)iV|. The perturbative behavior discussed in Sec. [il, based on the MS scheme, is 
qualitatively unchanged by the modification of the three-loop coefficient. For Nf = 12, the 
three-loop SF coupling has a fixed point at g 2 ~ 5.18, compared with g 2 w 5.47 at three 
loops in the MS scheme. At three loops in the SF scheme, the value of the critical exponent 
Eq. g is 7 = 0.286. 

The four-loop coefficient in the SF scheme has not yet been computed. But the fact 
that in the MS scheme the four-loop correction shifts the fixed point by less than 10% from 
its three-loop value suggests that the same may be true in the SF scheme. This and the 
relative smallness of the expected loop expansion parameter 0(^ 2 (L)/(47r 2 )) at the fixed 
point indicates that perturbation theory could be reliable to describe infrared behavior for 
Nf = 12, and that the infrared fixed point might truly exist. For Nf = 8, since the universal 
one- and two-loop coefficients are both positive, there is no reliable, perturbative evidence 
for the existence of an infrared fixed point. As observed in Sec. [TTj a nonperturbative study 
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is essential. 

In perturbation theory, the SF running coupling behaves just like the running coupling 
defined in other, more familiar ways. Its behavior is identical through two loops, and then 
qualitatively similar at three loops and beyond. Other definitions of the running coupling 
are based on Green functions of local operators or quantities such as the Wilson loop, while 
the Schrodinger functional and the associated running coupling dependence on a background 
field act across the entire lattice. This definition of the running coupling is nonperturbative, 
but its relation to other nonperturbative definitions in the strong-coupling regime is not yet 
clear. Nevertheless, the SF running coupling is adequate for our purposes: to distinguish 
between conformal and confining behavior in the infrared. 



IV. LATTICE SIMULATIONS 
A. Setup and Procedure 

To measure the running coupling on the lattice, we generate an ensemble of gauge config- 
urations distributed with the appropriate weighting by the Euclidean action. The running 
coupling is then computed as an average over this ensemble. Simulations are performed 
using the MILC code with some customization. Evolution of the gauge configurations 
is accomplished using the hybrid molecular dynamics (HMD) approach, with the fermionic 
contribution included using the R algorithm [25]. Trajectories are taken to be of unit length, 
while the step size At of the MD integrator is varied. The R algorithm is known to introduce 
errors of 0((Ar) 2 ) into observables; we discuss this effect along with other sources of error 
below. 

Sets of gauge configurations are generated at a fixed box size L/a and fixed bare coupling 
(3. For each ([3, L/a), two independent ensembles are created at T/a — L/a ± 1, and then 



averaged together as in Eq. (20). The data are collected over a wide range of (3 values and 
for 6 < L/a < 20, in order to capture the evolution of g 2 (L) over a large range of scales. It 
is important to note that in the range of (3 values employed, for both Nf = 8 and Nf = 12, 
there is no evidence for a bulk phase transition. We have explored this issue by examining 
the plaquette time series within this range and at lower values of (3. At lower values, we have 
indeed found evidence for a bulk phase transition. These lower values are, however, well 
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separated from the minimum (3 shown in our data tables and used in our analysis. There is 
no such evidence in the range of /3 values employed here. 

Simulations were also performed at L/a = 4, but these values were not used in our 
analysis. Examination reveals large lattice-artifact corrections in the L/a = 4 data. Their 
presence is not unexpected, particularly on the smaller 4 3 x 3 lattices; with only a single 
lattice site between the Dirichlet boundaries, the 0(a) boundary operators appearing in 
Eq. ^ overlap. In addition, the 4 3 x 3 lattices fail to satisfy the conditions of the "stability 
theorem" of Luscher et al, meaning that the background chromoelectric field being expanded 
around may not be an absolute minimum of the action [2]. This also precludes the use of 
data from 4 3 x 5 lattices, since although they satisfy the stability criterion, without the 
4 3 x 3 lattices we cannot use the averaging procedure described in Sec. 

Since updating of the gauge fields is accomplished locally, while the running coupling is 
simulated on the scale of the box size L, a large number of updates is required to generate 
statistically independent values of g 2 {L). To remove statistical bias from our results, we 
collect a large number of gauge configurations at each point in parameter space, ranging 
from 20000 to 160000 MD trajectories with a greater number required for measurements at 
stronger coupling. These run lengths are established based on our analysis of autocorrela- 
tions, discussed in Appendix [A| 

B. Step scaling 

We are interested in mapping out the behavior of the running coupling over a large range 
of scales, from the ultraviolet to the infrared. Often, a lattice simulation (i.e. a set of gauge 
configurations, generated using a fixed set of parameters) is focused on measuring quantities 
at distance scales lying between the lattice spacing a and the box size L. Our use of the 
Schrodinger functional instead places the observable g 2 (L) at the scale L, eliminating the 
latter restriction. However, the range of scales L over which we can measure the coupling 
strength with fixed lattice spacing a before the computational expense becomes prohibitive 
is still rather limited. To achieve our goal, therefore, we must measure the coupling using a 
wide range of a values, and then match these measurements together. To accomplish this, 
we use a procedure known as step scaling (26J, |27] . 

Step scaling provides a systematic way to combine multiple lattice measurements of the 
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running coupling g 2 (L) into a single measurement of the continuum evolution of the coupling 
as the scale changes from L — > sL, where s is a scaling factor called the step size. In a 
continuum setting, the evolution is described by the "step-scaling function," 



which can be thought of as a discrete version of the usual continuum evolution described 
by the beta function. In a lattice calculation, the extracted step-scaling function will be a 
function also of a/L, which we must extrapolate to the continuum: 



Step scaling is generically implemented on the lattice as follows. First, an initial value 
u = g 2 (L) is chosen. Several ensembles with different values of a/L are then generated, with 
f3 tuned so that the coupling measured on each is equal to the chosen value, g 2 (L) = u. 
A second ensemble is generated at each [3, but with L — > sL. The value of the coupling 
measured on this larger lattice is exactly E(s, u, a/L). An extrapolation a/L — > can then 
recover the continuum value a(s,u). Taking a(s,u) to be the new starting value, we can 
then iterate this procedure until we have sampled g 2 (L) over a large range of L values. In 
practice we take s = 2. 

There is a natural caveat on the step-scaling procedure. In the limit a/L — > with g 2 (L) 
fixed, g 2 (a/L) depends on the short-distance behavior of the theory, and it is important that 
it remains bounded so that it does not trigger a bulk phase transition. If asymptotic freedom 
governs the short distance behavior, this is automatic since g 2 (a/L) — > l/log(L/a). While 
this is our principal focus, the existence of an infrared fixed point for the Nf = 12 theory will 
lead us to consider also values of g 2 (L) lying above the fixed point. Then g 2 (a/L) increases 
as a — > 0, with no evidence from our simulations that it remains bounded and therefore that 
the continuum limit exists. Nevertheless, one can consider small values of a/L providing 
that (7q(o/L) remains small enough not to trigger a bulk phase transition. We return to this 
point in our discussion of the Nf = 12 simulation data. 

C. Interpolating Functions 

Carrying out the above procedure directly can be expensive in computational power since 
each tuning of (3 may require several attempts. The procedure also severely limits the rate 



a(s,f(L)) = f(sL), 



(22) 




(23) 
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at which computations may be performed, since each simulation must be finished and the 
value of a(s,u) extracted before the next iteration. We instead measure g 2 (L) for a limited 
set of values for (3 and L/a, and then generate an interpolating function. This function is 
then used to tune (3 as described above, and renders the cost of extracting a step-scaling 
function independent of the number of steps taken. 

For any value of L/a in our range, g~ 2 {(3,L) is a monotonically decreasing function of 
(3 = Q/g 2 . One procedure is simply to interpolate linearly between the available (3 values 
for each L/a. This works reasonably well in the perturbative region, reproducing the cor- 
rect continuum perturbative running once the step-scaling procedure is carried out. For 
stronger coupling, however, linear interpolation leads to some anomalies due to statistical 
fluctuations. A better procedure is to use a smooth interpolating function fit to the data. 

For our results reported in Ref . [1] , we fit g 2 ((3, L) to a single function based on a truncated 
Laurent series with poles at small, unphysical values of (3, well below the simulation range. 
Here we employ a set of interpolating functions, one for each L/a, focused directly on the 
lattice observable l/g 2 {(3, L/a). Motivated by the fact that in lattice perturbation theory 
this quantity takes the form 



we use a fit to g 2 (/3,L/a) at each L/a as a function of (3, with n-th order polynomial 
dependence on g 2 = 6/(3: 



The order n of the polynomial is varied with L/a in order to achieve the optimal x 2 per 
degree of freedom when fitting to our data. The values of the parameters with associated 
errors, determined by fits to the simulation data for both Nf = 8 and Nf = 12, are recorded 
in Appendix [Bj 

This function is used for interpolation within the measured range, as a basis for the 
step-scaling procedure. It is based on empirical observation of the g 2 dependence of our 
observable, and is not meant to imply that perturbation theory is applicable to our nonper- 
turbative, strong-coupling results. More elaborate interpolating functions could be used, in 
particular, modeling explicitly the L/a dependence or including nonanalytic terms in g 2 , but 




(24) 




(25) 



14 



such functional forms do not significantly alter the fit quality or the results of step scaling 
based on the collected data set. 



D. Statistical and systematic errors 

We account for numerous sources of statistical and systematic error in our analysis. A de- 
tailed discussion of the estimation and/or elimination of these errors is given in Appendix |Aj 
We conclude that potential systematic errors in our procedure are small compared to current 
statistical errors. Our final results for continuum running are therefore shown with only a 
statistical-error band. 

We note that this conclusion differs from that of Ref. [I]. In that reference, statistical 
errors were estimated in a less sophisticated way, in particular, ignoring the accumulation 
of error over repeated step-scaling steps. This led to very small statistical error bars. In 
contrast, the systematic error as determined by uncertainty in the correct form of the contin- 
uum extrapolation was large, due primarily to the inclusion of values of g 2 {L) computed on 



L/a = 4 volumes, which we now discard for reasons discussed in Sec. IV A Thus, with the 



more extensive analysis described here, statistical errors dominate rather than systematic. 



V. RESULTS 



Nf 



The simulation data for g 2 (L) as a function of f3 and L/a are displayed in Table [IJ Each 



data point is the average given by Eq. (20), with the statistical error in parentheses. The 
table ranges from (3 = 4.55-192 and L/a = 6-16. The lower limit is chosen to insure 
that the lattice coupling is weak enough so as not to induce a bulk phase transition. The 
upper limit is taken to be large so that we can check the agreement of our simulations with 
perturbation theory when the coupling is very weak. The final results reported here depend 
sensitively only on simulations below /? = 10. The data becomes more sparse with increasing 
L/a, reflecting the growing computational time involved. In particular, only a very limited 
amount of L/a = 20 data, at very weak coupling, is available at Nf = 8, so we exclude 
these points from our analysis. The L/a = 10 data is thus also excluded, since it cannot be 
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used in step scaling at s = 2 without the L/a = 20 points. The listed values of g 2 (L) are 
perturbative (g 2 (L)/4n <C 1) throughout much of the table, except for small f3. 

In order to carry out the step-scaling procedure, we employ the interpolating function of 



Eq. (25). The resulting best-fit mean values and errors for the parameters at each L/a are 
shown in Table [TTJ More details, including full covariance matrices, will be made available 
in the AIP's Electronic Physics Auxiliary Publication Service [28]. In Fig. [TJ data points 
are shown together with the interpolating functions for g 2 (L) as a function of /3, for each of 
L/a = 6, 8, 12,16. 

We implement the step-scaling procedure and extrapolation to the continuum as described 
in Sec. IV Figure [2] shows a typical continuum extrapolation from our 8-flavor data. The 
points shown represent steps from L/a = 6 — > 12 and 8 — > 16. Constant extrapolation 
(a weighted average of the two points) is used since the lattice-artifact contributions to 
E(2, u,a/L) are small compared to the statistical errors. We have estimated the systematic 
error in this procedure and found that it is small compared to the statistical error; details 
of this analysis are provided in Appendix [Aj 

Our results for the continuum running of g 2 {L) are shown in Fig. [3] We take L$ to be the 
scale at which g 2 (Lo) = 1.6, a relatively weak value. The points are shown for values of L/Lq 
increasing by factors of 2. The (statistical) errors are derived as described in Appendix |Aj 
For comparison, the perturbative running of g 2 (L) at two loops and three loops is shown 
up through g 2 (L) « 10 where perturbation theory is no longer expected to be accurate. 
The results show that the coupling evolves according to perturbation theory up through 
g 2 {L) « 4, and then increases more rapidly, reaching values that clearly exceed typical 
estimates of the strength required to trigger spontaneous chiral symmetry breaking [29J. 
The dynamical fermion mass is of order of the corresponding 1/L, and since the coupling is 
strong here, the theory will confine at roughly this distance scale. There is no evidence for 
an infrared fixed point or even an inflection point in the behavior of g 2 {L). 



B. N f 



12 



The simulation data for g 2 (L) as a function of (3 and L/a are displayed in Table III As 
with Nf = 8, each data point is the average given by Eq. (20), with the estimated error in 



parentheses. The table ranges from j3 = 4.2-192 and L/a = 6-20. The lower limit on j3 
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g 2 (L) 




FIG. 1: Measured values g 2 (L) versus for Nf = 8. The interpolating curves shown represent 



the best fit to the data, using the functional form Eq. (25). The errors are statistical, derived as 
discussed in Appendix |A} 

insures that the lattice coupling is weak enough so as not to induce a bulk phase transition. 
As in the Nf = 8 case, the upper limit is taken to be large in order to explore agreement 
with perturbation theory, but data above (3 = 10 do not have significant influence on our 
analysis. L/a = 20 data are included here and not in the Nf = 8 case because of concerns 
about the magnitude of the lattice artifact corrections, compared to the continuum running. 
In the end, artifact corrections were found to be small compared to our statistical error. 
Data become more sparse with increasing L/a, reflecting the growing computational cost 



involved. The interpolating functional form Eq. (25) is again employed, and the resulting 
best-fit mean values and errors of the parameters at each L/a are shown in Table IV The 
full covariance matrix will be made available in the AIP's Electronic Physics Auxiliary 
Publication Service [28J. In Fig. |4j data points are shown for g 2 (L) as a function of /3, 
together with the interpolating functions for each of L/a = 6, 8, 10, 12, 16, 20. 

The data and the interpolating curves already suggest the existence of an infrared fixed 
point for Nf = 12. For small /3, the general trend is that g 2 {L) decreases with increasing L. 
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FIG. 2: Step-scaling function E(2, u, a/L) at various u, for each of the two steps L/a = 6 — > 12 and 
8 — > 16 used in the TVy = 8 analysis. Note that H(2,u,a/L) > it in each case, with the difference 
increasing as u increases. 

This behavior and the fact that for larger (3, g 2 {L) increases with increasing L, is reflected 
in a crossover behavior in the interpolating curves. We first implement the step-scaling 
procedure choosing an initial u = g 2 (L) well below a possible fixed-point value so that a 



continuum limit is guaranteed to exist, as discussed in Sec. |IVB 

A constant continuum extrapolation (a weighted average of the available values of 
S(2,u,a/L)) is again employed for each u. Now, since we have data at L — 20, the ex- 
trapolation is a weighted average of three points corresponding to the steps 6 — >■ 12, 8 — ^ 16, 
and 10 — > 20. Examples of such a continuum extrapolation are shown in Fig. |5j The 
systematic error is again estimated to be small compared to the statistical error. 

Our results for the continuum running of g 2 {L) from below a possible infrared fixed point 
are shown in Fig. [6j L is again taken to be the scale at which g 2 (L ) = 1.6, a relatively 
weak value. The points are shown for for values of L/Lq increasing by factors of 2. The 
(statistical) errors are derived as described in Appendix [Aj For reference, the two-loop and 
three-loop perturbative curves for g 2 (L) are also shown in Fig. [6j From the figure, we 
conclude that the infrared behavior is indeed governed by a fixed point whose value lies 
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FIG. 3: Continuum running for Nf = 8. Purple points are derived by step-scaling using the 
constant continuum-extrapolation of Fig. [2j The error bars shown are purely statistical, and are 
derived as described in Appendix [A] Two-loop and three- loop perturbation theory curves are shown 
for comparison. 

within the statistical error band. Because of the underlying use of an interpolating function, 
the error bars of adjacent points in Fig. [6] are highly correlated. As the running coupling 
approaches the infrared fixed point, this correlation approaches 100%, so that the error bars 
asymptotically approach a stable value as the number of steps is taken to infinity. The range 
of possible values of the fixed point from our simulations is consistent with the three-loop 
perturbative value in the SF scheme. It is well below estimates [2H] of the strength required 
to trigger spontaneous chiral symmetry breaking and confinement. 

The infrared fixed point also governs the L — > oo behavior starting from values of g 2 {L) 



above the fixed point. As discussed in Sec. |IVB[ the continuum limit is then no longer 
guaranteed to exist and the step-scaling procedure cannot be naively applied. Instead, one 
can restrict the discussion to finite but small values of a/L, small enough to minimize lattice 
artifacts but large enough so that for g 2 (L) near but above the fixed point, g^a/L) remains 
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FIG. 4: Measured values g 2 (L) versus (3, Nf = 12. The interpolating curves shown represent the 



best fit to the data, using the functional form of Eq. (25). 
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FIG. 5: Step-scaling function E(2,u, a/L) at various u, for each of the three steps L/a = 6 — » 12, 
8 — > 16, 10 — > 20 used in the iVj = 12 analysis. Note that S(2, u, a/L) — > u as the starting coupling 
u approaches the fixed point value. 
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FIG. 6: Continuum running for Nf = 12. Results shown for running from below the infrared fixed 
point (purple triangles) are based on g 2 {Lo) = 1.6. Also shown is continuum backwards running 
from above the fixed point (light blue squares), based on g 2 (Lo) = 9.0. Error bars are again purely 
statistical, although strongly correlated due to the underlying interpolating functions. Two-loop 
and three-loop perturbation theory curves are shown for comparison. 

small enough not to trigger a bulk phase transition. Since we use a constant extrapolation, 
this procedure can be taken to define, within our errors, a g 2 {L) at a small but finite a/L. 
The step-scaling procedure then leads to the continuum running from above to the fixed 
point, also shown in Fig. [6j The statistical-error band is derived as in the approach from 
below. 

Finally we note that the exponent 7 governing the approach to the infrared fixed point 
in the SF scheme can also be extracted from the simulation data. Taking the log of Eq. Q, 
we see that the quantity log \g* — g 2 (L)} should have a linear dependence on L with slope 
—7 near the fixed point. Computing this quantity from our data, running from either above 
or below the fixed point, we find 7 = 0.13 ± 0.03, somewhat smaller than the three-loop SF 
perturbative estimate of 0.286. 
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C. Comparison with Other Lattice Work 



1. Schrodinger functional studies 

Lattice simulations for the SU(3) Schrodinger functional running coupling have been per- 
formed for Nf = 16 [30], for the quenched theory [31], and for Nf = 2 [4J. For Nf = 16 
[30] . the perturbative infrared fixed point is very weak. In this case, the simulations were 
done for values of the lattice coupling in the weak-coupling (chirally symmetric and decon- 
fined) phase but leading to values of g 2 {L) well above the perturbative fixed point. Evidence 
was presented that g 2 (L) decreases with increasing L, consistent with the approach to the 
fixed point from above as expected with a continuum infrared fixed point. A continuum 
extrapolation via the step-scaling procedure was, however, not implemented. 

For both the quenched theory [32] and Nf = 2 [1], the step-scaling procedure was im- 
plemented and a continuum running coupling was extracted. In each case, starting with a 
g 2 (L) well into the perturbative regime, the coupling grows through large, nonperturbative 
values. And in each case, the growth is more rapid than for Nf = 8, as shown in Fig. [3j For 
the quenched theory, g 2 (L) was argued to grow exponentially at large L, consistent with the 
leading order prediction from the strong-coupling expansion [32J. 

2. Other multifermion studies 

Lattice simulations of SU(3) gauge theories with multiple fermions in the fundamental 
representation began more than 20 years ago. Brown et al. [33] examined the Nf = 8 
case with staggered fermions, providing evidence that the theory confines, but remaining 
inconclusive due to finite volume effects. Damgaard et al. [M] examined the staggered 
Nf = 16 case, noting that even with the very weak infrared fixed point present in this 
theory, a bulk chiral transition sets in at sufficiently strong lattice coupling. In a 2001 
Columbia PhD thesis [6], Sui studied QCD for Nf = 2 and Nf = 4 staggered fermions, 
observing stronger finite-lattice-size effects in the latter case. The work of Iwasaki et al. [7] 
is perhaps most directly relevant to the results reported here and in Ref . [TJ . Through a focus 
on the strong lattice-coupling phase using Wilson fermions, they concluded that 6 < Nf < 7, 
in disagreement with our results. 

Interest in multifermion studies has grown considerably in the past few months. Deuze- 
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man et al. [35] have examined chiral symmetry breaking for the Nf = 8 case using staggered 
fermions, concluding that the lower end of the conformal window is indeed above Nf = 8. 
Jin and Mawhinney [35] have come to the same conclusion through a study of the chiral 
condensate and the heavy quark potential. Fodor et al. [37] have begun a multifermion 
simulation using staggered fermions, while Bilgici et al. have developed a new approach to 
running coupling measurement with an eye towards eventual multifermion measurements 



VI. SUMMARY AND DISCUSSION 



We have concluded from lattice simulations of the Schrodinger functional running cou- 
pling that for an SU(3) gauge theory with Nf Dirac fermions in the fundamental representa- 
tion, the value Nf = 8 lies outside the conformal window, leading to confinement and chiral 
symmetry breaking, while Nf = 12 lies within the conformal window, governed by an in- 
frared fixed point. We have bounded the fixed point value as shown in Fig. [6] and estimated 
the exponent 7 describing the approach to the fixed point Eq. This is, as far as we 
know, the first nonpertubative evidence for the existence of infrared conformal behavior in 
a nonsupersymmetric gauge theory. These results confirm and refine the analysis of Ref. pQ . 

The Nf = 8 and Nf = 12 results imply that the lower end of the conformal window, 
Nf, lies in the range 8 < Nj < 12. This conclusion, in disagreement with Ref. [7j, is 
reached employing the Schrodinger functional (SF) running coupling, g 2 (L), defined at the 
box boundary L with a set of special boundary conditions. This coupling is a gauge invariant 
quantity, valid for any coupling strength and running in accordance with perturbation theory 
at short distances. 

For Nf=8, we have simulated g 2 (L) up through values that exceed typical estimates of 
the coupling strength required to trigger dynamical chiral symmetry breaking [29], with no 
evidence for an infrared fixed point or even an inflection point. For Nf— 12, our observed 
infrared fixed point is rather weak, agreeing within the estimated errors with the three-loop 
fixed point in the SF scheme, and well below typical estimates of the coupling strength 
required to trigger dynamical chiral symmetry breaking [29J. 

Whether perturbation theory can be used reliably to reproduce the behavior in the vicin- 
ity of the Nf = 12 fixed point remains to be seen. The three-loop value of the fixed point 
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is substantially different from the two-loop value. On the other hand, in the MS scheme 
where the four-loop beta function has been computed, the four-loop fixed point is shifted by 
only a small amount from the three-loop value. The relative weakness of this fixed point, 
together with the fact that Nf cannot be much smaller, raises the question of whether the 
theory remains perturbative throughout the conformal window as suggested by Gardi and 
Grunberg [T2]. If this is the case, the behavior in the neighborhood of Nf would be rather 
different from the supersymmetric SU(iV) gauge theory [15] . In particular, it is not clear 
whether there would be a useful, effective low energy description of the infrared behavior. 

It is important to confirm our results by employing other definitions of the running 
coupling, for example, based on the Wilson loop and static potential [Sj, and by examining 
scheme-independent quantities. Most notably, spontaneously chiral symmetry breaking as a 
function of Nf must be studied through a zero-temperature lattice simulation of the chiral 
condensate [36]. Simulations of g 2 (L) for other values of Nf, in particular Nf— 10, are crucial 
to determine more accurately the lower end of the conformal window and to study the phase 
transition as a function of Nf. All of these analyses should be extended to other gauge 
groups and other representation assignments for the fermions [381 ESS KB EJ H2J H31 EU l4"5l 

SB snug. 

The phenomenological relevance of these studies remains to be seen. One possibility is 
that a theory with infrared conformal symmetry could describe some new sector, coupled 
to the standard model through gauge-singlet operators [19]. Another possibility, much dis- 
cussed in the literature, is that a theory with Nf outside but near the conformal window 
(< NJ) could describe electroweak breaking and provide the basis for walking technicolor 
[5TH loTl |5"2"] . In this class of theories, as Nf — > Nf from below, a hierarchy emerges be- 
tween the electroweak scale and the larger mass scale where the gauge coupling becomes 
strong. This could be signaled by the appearance of a plateau of finite extent in g 2 (L), 
and by the development of a hierarchy between the chiral condensate and the electroweak 
scale. It is also important to explore the particle spectrum in this limit and to compute 
the electroweak precision parameters, in particular the S parameter |53j. These studies are 
currently underway [9]. 
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APPENDIX A: STATISTICAL AND SYSTEMATIC ERROR IN STEP-SCALING 
ANALYSIS 

1. Numerical-simulation error 

The hybrid molecular dynamics (HMD) method involves the solution of classical equa- 
tions of motion, requiring a numerical integration at a finite step size At. This introduces 
systematic numerical errors of 0((Ar) 2 ) into all observables, including the running coupling. 
Removal of this finite step-size error from a given measurement of g 2 (j3,L) can be accom- 
plished by simulating at multiple values of At and performing a quadratic extrapolation to 
zero. At relatively weak values of the bare coupling, the step-size error is observed to reduce 
the measured values of g 2 . Furthermore, the magnitude of the shift increases with the box 
size L/a. Although statistically significant, the effect on step-scaling results is negligible for 
L/a = 8 and below, but becomes significant at larger box sizes. At stronger values of the 
bare coupling, no systematic shift due to step-size error can be resolved within our statistical 
error. Wherever the effect or step-size error is observed to be significant, only extrapolated 
values of g 2 are used in our analysis, effectively removing this source of systematic error. 

Although our goal is to generate a set of statistically independent gauge configurations, 
in practice configurations which are separated by only a small number of updates (by a small 
MD time) can be correlated with each other. The existence of these "autocorrelations" can 
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lead to underestimation of statistical errors if classical estimators are used. Even with appro- 
priate statistical methods, the time series must extend to several times the autocorrelation 
time of the observable in order to obtain an unbiased measurement. Statistical estimates 
of the integrated autocorrelation time show a clear trend towards longer autocorrelations 
as either the box size L/a or the strength of the lattice coupling g\ is increased, with the 
longest estimated time on the order of 1000 trajectories. 

In addition, autocorrelation times in certain time series are, in effect, enhanced by the 
observed phenomenon of "excursions." This has been noted in prior studies of the SF running 
coupling [221 The coupling is seen to jump to a new equilibrium value, which can remain 
stable for up to several thousand trajectories. This can be interpreted as a tunneling of the 
system from near the original minimum of the action, determined by the imposed background 
field, into other metastable minima. In the presence of these excursions, autocorrelations 
in the observable are introduced on the scales of the average duration and period of the 
tunneling events. 

Although the excursions themselves have a physical interpretation, the correlations that 
they induce in the data are artifacts of the procedure used to generate gauge configurations, 
and the associated time scales are dependent on the choice of such a procedure. For example, 
an update algorithm based on the selection of completely random gauge configurations would 
have no autocorrelations by design, but could still show evidence of tunneling into secondary 
minima in the form of non-Gaussianity in a histogram of measured coupling strengths. 

The excursions are empirically seen to occur always to an equilibrium at stronger coupling 
g 2 {L), and become more frequent as the bare coupling strength is increased from weak to 
stronger values. In particular we observe tunneling events in the running average of the 
time series for 5.8 < (3 < 4.7, at both eight and twelve flavors, with these events becoming 
impossible to isolate from statistical fluctuations at stronger coupling. Some representative 
plots showing this effect are shown in Fig. [7} The contrast between the two time series is 
sharp, with tunneling events clearly evident only in the stronger-coupling time series. 

We therefore choose the target number of trajectories for a particular measurement of 
g 2 {(3,L) as follows: 20000 trajectories for (3 > 8.0; 40000 for 5.0 > (3 > 8.0; and 80000 for 
(3 < 5.0. These values exceed all above estimates of the autocorrelation time, and allow 
sampling of multiple excursions where such events are observed. 

Estimation of statistical error, with full propagation of errors including continuum- 
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FIG. 7: Representative time series taken from our 12-flavor simulations. The estimated mean 
values with statistical error for each time series are depicted in the dotted lines and the points with 
error bars to the right of the plot. The time-series plot is a running average over a window of 800 
points. Both data sets were gathered on volumes of 16 3 x 15, and represent relatively weak coupling 
{(3 = 6.5, light blue) and somewhat stronger coupling {(3 = 4.7, dark red). The difference in length 
between the two time-series reflects our choice of target numbers of trajectories, as discussed in 
the text. 



extrapolation error through all step-scaling steps, is performed using the bootstrap method. 
The raw data are first reduced to uncorrelated blocks. Two thousand bootstrap replications 
of the data set are generated, and quantities of interest are computed as statistics on the 
bootstrap data. Two-sided errors are shown in all cases, representing la confidence inter- 
vals on the mean, computed using the bias-corrected and accelerated confidence interval 
estimation method [56J. Measured values for g 2 {L) with estimated two-sided error bars are 
included in Tables [T] and III Although we quote a single value for each (/3, L/a) in the data 
tables, we emphasize that our analysis is carried through on the complete sets of bootstrap 
replicated data from which these mean values are derived. 
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2. Interpolating- function error 



The choice of a particular interpolating functional form is a potential source of systematic 
error, particularly if it yields a poor fit to the data, reflected by a value of x 2 P er degree of 



freedom (DOF) significantly larger than 1. Using our interpolating functional form Eq. (25), 
we find an excellent fit to the simulation data; at Nf = 8, we find an overall x 2 /DOF= 
1.62 ± 0.30 with 107 degrees of freedom, while for the Nf = 12 data the fit yields x 2 /DOF= 
1.47 ± 0.26 with 171 degrees of freedom. Errors on x 2 are estimated using the jackknife 
method. Note that these are global x 2 values, representing a sum over contributions from all 
L/a fits, divided by the total number of degrees of freedom. The relatively low probability 
of these values in the \ 2 distribution indicates that our error bars are likely somewhat 
underestimated. Possible sources of this effect are detailed in the previous subsection. 

We have further attempted to search for systematic error due to the interpolating function 
by trying to fit a variety of other functional forms. One approach is to attempt to improve 



upon the form in Eq. (25) by the addition of extra terms. We have considered the addition 
of terms nonanalytic in g%, and the addition of terms and constraints which reproduce 
one-loop perturbation theory in the limit Qq — *■ 0. In each case, the additional terms do 
not significantly improve x 2 /DOF, and the results of the analysis based on such fits are 



indistinguishable from those based on Eq.(25). 

Another possibility is to use an altogether different functional form, such as the Laurent 
series expansion of Ref. |TJ . These forms can also be extended with nonanalytic terms and 
perturbative constraints. For Nf = 12, large systematic shifts are seen in the continuum 
running curves based on such fits. This systematic effect is universally associated with a 
significantly higher x 2 /DOF, indicating that the interpolating function does not accurately 
reflect the underlying measurements. 

We stress that since our fit functions are used only for interpolation, not extrapolation, 
any two fits to the coupling measurements which yield comparable and acceptably small 
X 2 /DOF will give indistinguishable results for the continuum running. 

We conclude that the systematic error associated with the selection of the final form 



Eq. (25) is negligible, given the values of x 2 /F>OF quoted above. 
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3. Continuum-extrapolation error 

Since each step in the step-scaling procedure involves a continuum extrapolation of the 
step-scaling function H(2,u,a/L), the choice of functional form in the extrapolation is an- 
other potential source of systematic error. 

What is the expected behavior of £(2,w, a/L) as a function of a/L? Staggered fermions 
with the Wilson gauge action inherently lead to discretization effects of 0(a 2 ). However, the 
presence of the Dirichlet boundaries here leads to operators which contribute 0(a) lattice 
artifacts, as shown in Eq. Q. We include a counterterm for these operators, with its value 
determined by one-loop perturbation theory. With this counterterm, we expect that the 
0(a) terms are small, even more so in E(2,u, a/L) where much of the a/L dependence is 
removed. 

By comparing our data to perturbative running at fixed, relatively weak lattice coupling, 
we find that in fact all lattice-artifact corrections are negligible compared to our statis- 
tical errors for Lja > 8. Independent of perturbation theory, we find that the best fit 
(smallest x 2 /DOF) for the continuum extrapolation at Nf = 12 is given by a constant, a/L- 
independent extrapolation, yielding <r(2,u) as a weighted average of S(2,w, a/L) over the 
available range of a/L. Based on this experience, we use constant extrapolation at Nf = 8 
as well. Lattice artifacts play a less important role here because of the stronger continuum 
running (which also provides justification for not including data at L/a = 20). The absence 
of Lja = 20 data at Nf = 8 means that there are only two available steps (6 — > 12 and 
8 — > 16). Thus the constant a/L extrapolation is the only constrained fit available. The 
errors for each T,(2,u,a/L) lead to statistical errors in a(2,u), which is represented by a 
statistical-error band for the final continuum-running curve. 

Statistical error in the continuum extrapolation is computed by the bootstrap method; 
the extrapolation is performed independently on each bootstrap ensemble, leading to a 
distribution of values for a(2,u), which can then be used to estimate a mean value and 
two-sided confidence interval. The result of the application of n steps, a(2 n ,u), is likewise 
computed within each bootstrap ensemble to obtain a full distribution. 
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APPENDIX B: RUNNING COUPLING DATA 



Our measurements of the running coupling g 2 (L) are presented in Tables [T] and III below 



Two-sided error bars are estimated using the BCa method as described in Appendix [Aj 
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TABLE I: Measurements of g 2 (/3, L/a), N f = 8. 
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9.5(+l"4) 

V — l.z / 


4.900 


4.761 

V OZ / 


5.18(1^) 

V If/ 








5.000 


4-204 (+11) 

\ — DO / 


4.68(^|) 

V — 10 / 


4-96(+lD 

V — Iz / 


5.64(+34) 


6.45 (+5J) 

v — 4y / 


5.100 


3.788(li) 

V O.J / 


4.2(±g-i) 

V U. -L / 








5.200 


3.382(H) 


3.674(1^3) 








5.300 


3.087 (+11) 


3.31l(+35) 








5.400 


2.972 (+^) 

v zy / 


3.145 (If*) 








5.500 


2. 724 (+22) 


2.965 (+1) 

V Oo / 


3.122(1S) 


3.380 (+90) 


3.374(+il) 

v — OO / 


5.600 


2.603(+ji) 


2.795 (+43) 

\ O-L / 








5.700 


2.4424 ( +1 ° ) 


2.590(+}3) 

V 10 / 








5.800 


2.3340(1|) 


2.491(1^) 

V J-O / 








5.830 


2.286(+ 9 1 ) 

V 11/ 


2.456(111) 

V 10 / 


2.535(+^) 

V IO / 


2.647(+|) 

V ZO / 


2.842 M\) 


5.900 


2.2246(1^) 


2-340(111) 








6.000 


2.1374(H) 


2.264(1{ 2 ) 








6.100 


2.0578(1^3) 


2.1622(1^) 








6.200 


1.9778(1^) 


2.0687(l|t) 








6.300 


1.9127(1|) 


2.0008(H) 








6.400 


1.8430(H) 


1.9230(+Jj>) 
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?{L) 


L/a 





6 


8 


10 


12 


16 


6.500 


1.7869 (+50) 


1.8601( + S) 








6.590 


1.7328 ( + f,) 

\ —61 ) 










6.600 


1.7232 (+39) 

V — oo / 


1.7996 (+|4) 

V — 4D / 








6.700 


1.6665 ( + f 4 ) 

V — 44 / 


1.7347( + jjl) 

V — DO / 








6.800 


1.6253 ( + ff.) 

\ —zo / 


1.6938(+S) 

V — OO / 


1.7359 

\ Do / 






6.883 


1.5865 (+^) 

v — ZL / 






1.7143( 


f82\ 

— to / 


6.900 


1.5776 (+3?) 

V — ol / 


1.6260( + §}) 

V —01 / 








7.000 


1.5323(+S2) 

\ — OO / 


1.5845 ( + ^) 

V — oo / 






1.734(+9) 

V —11 / 


7.090 






1.585(+}2) 

V — 16 ) 






7.100 


1.4886 (+22) 

V -^O / 




1.5882 (+6°) 






7.115 






1.5787 






7.153 






1.5664(+|) 

v —oo / 






7.200 


1.4559 ( +22) 

V —27 ) 










7.300 


1.4133( + |Z) 










7.400 


1.3795 (+22) 

v — Z ( / 










7.500 


1.3436 (+||) 

V — Z6 / 


1.3880( + 5i) 

V — oZ / 








8.000 


1.2000(+^) 

v — iy / 


1.2302 

\ — OO / 








8.500 


1.0863 (tJI) 


1.1141(124) 








9.000 


0.99470 

\ — yo / 


1.0189( +2 H) 

V — Zif 








12.000 


0.65805 (+11) 


0.66857(+n) 




0.6834( 


tH) 0.7054(+68) 

O-L / \ — UO / 


12.800 


0.6036 ( +, 9 n ) 

v — 1U / 


0.61228(+2S) 

v — yy / 








13.710 


0.55193(+fZ) 


0.5618(+ft) 

V — 11 / 








14.770 


0.50383( + !o) 

V —48 / 


0.50818( + H) 

\ — o4 / 








16.000 


0.45600(1^) 


0.46117(H) 








17.450 


0.41039(1^) 


0.41453(1|^) 








19.200 


0.36557(1^) 


0.36869(1^) 








21.330 


0.32349(1^) 


0.32680(1^) 








24.000 


0.282902(1^) 0.28465(li|) 




0.28876( 


-78 ) 
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?{L) 


L/a 





6 8 


10 


12 


16 


27.430 
32.000 
38.400 
48.000 
64.000 
96.000 
192.000 


0.24331 (tig) 0.24490(+J6) 
0.2054O(ilf) 0.20632(1|) 
0.1686l(+}£) 0.16893(+i) 
0.13248( +! 9 ) 0.133007(1^) 
0.097919(1^) 0.098113(1^) 
0.064344(1^) 0.064433(1^) 
0.031717(11|) 0.031675(1||) 




0.1336l(li) 





TABLE II: Interpolation best-fit parameters, Nf = 8. 





L/a 


par am 


6 


8 


12 


16 


c l,L/a 


0.4632(99) 0.4932(59) 0.581(18) 


1.01(18) 


C2,L/a 


-0.14(12) 


-0.167(44) 


-0.235(40) 


-1.01(37) 


c 3,L/a 


1.13(59) 


0.76(11) 


0.245(23) 


0.60(19) 


c 4,L/a 


-3.2(1.3) 


-0.98(11) 








c 5,L/a 


4.7(1.6) 


0.441(37) 








c 6,L/a 


-3.43(92) 











c 7,L/a 


0.98(21) 











X 2 /dof 


1.78 


1.49 


1.33 


1.65 


N dof 


55 


43 


7 


2 
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TABLE III: Measurements of g 2 ((3, L/a), N f = 12. 



?{L) 


Lja 


a 


6 


8 


10 


12 


16 


20 


4.200 


14.84 

\ —ao / 


11.3l(+fJ) 

V — OO / 


i3.o(+i-n 

V — 1.0 / 


12.5(+11) 

V — 1.1 / 






4.210 


13.63( + 2Z) 

\ — ao / 


11.12(1$) 










4.220 


12.44( +21) 

V AO / 


ii.i(18I) 










4.230 


11.6(±g-l) 


8.9(±l-f) 










4.240 


n.io(i^) 

V If/ 


7.7(^1) 










4.250 


10.48 (+|) 

V AO / 


9.14(±i) 


9-32(1?°) 








4.270 


9.47(1}?) 


8-75(11?) 










4.300 


8.3(+°-l) 

v — U.l / 


7.42 (+30) 


7.97(+|) 

V — oo / 


7.49 (+41) 


8-38 ( + H) 

\ — yo / 


6-62 ( + H) 

V — 45 / 


4.350 


6.91(+1|) 

V —LO / 


6.84(+l) 


6.84 (+2?) 

V —Al / 


7.39 (+1) 

V —Oo / 






4.400 


6.140(1??) 


6.08(+30) 

V a\) j 


6-14( + |) 

V OO / 


5.95(H) 






4.450 


5.474 


5.37(li) 

V Z4 / 


6.00 (IP) 

V A 1 / 








4.500 


5-151(1^) 


5.05(111) 

V i-A / 


5.26(111) 


4 -97(l^) 






4.600 


4.248 (If ) 


4.46(111) 


4.42(H|) 


4.22(H7) 

V -LO / 


4-01(11) 




4.700 


3.822 (+11) 


3.746(111) 


3.669 (+i 7 2 4 °) 


3.81(11?) 

V A 1 / 


5.00(li) 




4.800 


3.458 (+33) 


3.549 (+1) 
v — yo / 


3.58(+l?) 

V — 11 / 








4.900 


3.061 (+ 4 ?) 


3.281 (+25) 


3.196(+fT) 

V Ol / 








5.000 


2.884(+i) 

V AA / 


2-912 (1|°) 

\ — AO / 


3.005 (+1) 

v — oO / 


3.023(+? 4 ) 

V 11/ 


3.28(+l) 

V AO / 




5.100 


2.733(+ 43 ,) 


2.852 (+11) 


2.811(+||) 

\ — AO / 








5.200 


2.549(+{3) 

V — 1 o / 


2.573(+i) 










5.300 


2.466(+i) 

V AO / 


2.438(li) 


2.528(1||) 

V — AO / 








5.400 


2.325(1|) 


2.337(11?) 










5.500 


2.1985 (t 1 ^ ) 


2.2346(+<+) 


2.273(111) 


2.27l(li|) 


2.360(131) 


2.311(1^) 


5.600 


2.0979(1^) 


2.148(113) 










5.700 


2.0139(1 4 J) 


2.066(111) 


2.1143(1??) 








5.800 


1.9461(1^) 


2.016(131) 










5.900 


1.8636(1|) 


1.8970(H) 


1-935(111) 
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L/a 


(3 


6 


8 


10 


12 


16 


20 


6.000 


1.8039(+^) 

V —44 / 


1.8218(+*°) 

V — o4 / 


1.879(+{}) 

V —11 / 


1.873( + }H 

V — lo / 


1.912( + "j) 

V — lo / 


1.922(+22) 

V —ZO } 


6.100 


1.7532( + 2Z) 

V —Zo / 


1.7862 (+^ 2 ) 

V — bl / 


1.813(+|}) 

V — 11 / 








6.200 


1.6909 (+48) 

V —44 / 


1.7400 (+^) 

V — oU / 










6.300 


1.6405(+?°) 

v — iy / 


1.6698 (+5?) 


1.6901 ( + f n ) 

V — oU / 








6.400 


1.5827(1^) 


1.6299 










6.500 


1.5459 (+^) 

V —24: ) 


1.5702 (+^ 2 ) 

V — ol / 


1.5967(+|i) 


1.600(+") 

V — 10 I 


1.614(+Ji) 

V — lo 1 




6.600 


1.5031 ( 

V —15 / 


1.5274 ( + IT) 

\ — 25 / 










6.800 


1.4253(+^) 

V — 14 / 


1.4509 ( +3 °) 

V — Ol / 










7.000 


1.3471 ( +2 .J) 

V —24 1 


1.407( +2 |) 

\ —26 1 


1.3985 ( + fo) 

\ — 45 / 


V — 18 ) 






7.200 


1.2932(+} 2 ) 

\ —IZ / 


1.3106 f^) 

V — ZO / 










7.400 


1.2354(+}{) 


1.2538 ( +2 1) 

V — zz / 










7.500 




1.2163( + «f) 

V — ol / 




1.2440 (+*!) 

V — OO / 


1.2465(+™) 

V —oo / 




7.600 


1.18292( 

V —95 / 


1.2000(1^) 










7.800 


V —so ) 


1.1516(+}2) 

V — 19 / 










8.000 


1.0880(+}5) 

V — 14 / 


1.1027(+£f) 


1.1197(+i£) 

V — oU / 


1.1332 (+ 120 ) 


1.1312(+|1) 

V — y4 / 


1.1415(+s°) 

V —o i / 


12.000 


0.62354( + ") 

V —DO / 


0.62833( + ^) 


0.63109 

v — yi / 








12.800 


0.57479 ( + f R ) 


0.5798 ( +, 9 n) 

V — 1U / 


0.58266 (+%) 








13.710 


0.52879 (+11) 

\ — ^y / 


0.53198(+S) 

V / u / 


0.53620 (+15) 








14.770 


0.48276 (+??) 


0.48656 

V Ol / 


0.48839 (+f 2) 








16.000 


0.43917( +3 i) 


0.44166(+5q) 

v — 4y / 


0.44383 (+11) 

V — o4 / 




0.4482 (+}?) 


0.4479 (+?§) 

V —ZD } 


17.450 


0.39659 

V — ol / 


0.39940(+^) 

v — oo / 


0.40084 (+60) 

v — oz j 








19.200 


0.35498 (+ 23 ,) 

\ —zo / 


0.35677( +3 ?) 

V — ol / 


0.35844(+$j) 

V — oo / 








21.330 


0.31519( + }|) 

V — lo / 


0.31732( +2 .£) 

V —29 / 


0.31837 ( 








24.000 


0.27629 ( 

■ y 1 ■ i / 


0.27779(11) 


0.27897 (±40) 


0.27852 






27.430 


0.23839(1^) 


0.23937(1^) 


0.24030 (iH) 


0.23987(1^|) 






32.000 


0.20157(1^) 


0.20222(1^) 


0.20323 (til) 


0.20257(1 22 ) 






38.400 


0.16603(1}|) 


0.16643(1^) 


0.16710(1|) 


0.16740(1^) 






48.000 


0.13127(1^) 


0.13124(1^) 


0.13178(1^) 


0.13146(1^) 
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f(L) L/a 

(3 6 8 10 12 16 20 

64.000 0.097222 (+^) 0.097401 ( l 1 ^ ) 0.097317(1^) 0.097257(1^) 
96.000 0.063998(1^) 0.064060(1^) 0.063967(1^) 0.064089(1^) 
192.000 0.031638 (1^) 0.031645(1^) 0.031677(1^) 0.031673(1^) 0.031675(1^) 



TABLE IV: Interpolation best-fit parameters, Nf = 12. 





L/a 


par am 


6 


8 


10 


12 


16 


20 


c l,L/a 


0.380(11) 


0.4092(66) 0.4269(97) 0.413(10) 


0.467(20) 


0.463(32) 


c 2,L/a 


-0.08(13) 


-0.192(46) 


-0.224(70) 


-0.167(84) 


-0.154(46) 


-0.111(70) 


c 3,L/a 


0.56(54) 


0.73(11) 


0.75(17) 


0.82(22) 


0.164(28) 


0.129(38) 


C±,L/a 


-1.2(1.1) 


-0.837(98) 


-0.81(17) 


-1.02(23) 








c 5,L/a 


1.6(1.1) 


0.342(31) 


0.319(54) 


0.417(80) 








c 6,L/a 


-1.10(57) 

















c 7,L/a 


0.32(12) 

















x 2 /dof 


1.60 


1.42 


1.30 


1.59 


1.61 


1.10 


Ndof 


54 


55 


36 


17 


7 


2 
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